Single molecule pulling with large time steps 
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Recently, we presented a generalisation of the Jarzynski non-equilibrium work theorem for phase 
space mappings. The formalism shows that one can determine free energy differences from approx- 
imate trajectories obtained from molecular dynamics simulations in which very large timesteps are 
used. In this work we test the method by simulating the force induced unfolding of a deca-alanine 
helix in vacuum. The excellent agreement between results obtained with a small, conservative time 
step of 0.5 fs and results obtained with a time step of 3.2 fs (i.e., close to the stability limit) indicates 
that the large time step approach is practical for such complex biomolecules. We further adapt the 
method of Hummer and Szabo for the simulation of single-molecule force spectroscopy experiments 
to the large time step method. While trajectories generated with large steps are approximate and 
may be unphysical — in the simulations presented here we observe a violation of the equipartition 
theorem — the computed free energies are exact in principle. In terms of efficiency, the optimum 
time step for the unfolding simulations lies in the range 1-3 fs. 
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I. INTRODUCTION 

Modern experimental techniques such as atomic force 
microscopy or optical tweezers offer the means to study 
the mechanical properties of single molecules such as pro- 
teins or nucleic acids providing striking insights into their 
structure, energetics and dynamics [J HI, H[- In such ex- 
periments, finely tuned forces are used to distort indi- 
vidual molecules and to track their response with high 
resolution. Although mechanical single molecule exper- 
iments can yield a wealth of useful information, their 
thermodynamic analysis in terms of binding constants 
or unfolding free energies is not straightforward. This 
complication originates from the fact that typically the 
perturbation acting on the system is time dependent, 
driving the system away from equilibrium. In such a 
non-equilibrium situation the free energy, or reversible 
work, cannot be simply calculated by integration of the 
average force along the pulling path. A solution is of- 
fered by Jarzynski's theorem which relates equilibrium 
free energies to the statistics of work carried out during 
non-equilibrium transformations 0, [f| . 

According to the Jarzynski equality, the free energy 
difference AF between two states corresponding to differ- 
ent values of an externally controlled parameter A can be 
obtained from an exponential average of the work W per- 
formed on the system by switching the parameter from 
its initial to its final value: 
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The average, denoted by the angular brackets, extends 
over all realisations of the switching process starting from 
a system initially in equilibrium with a heat bath at tem- 
perature T = l/k^T. In an experiment in which a single 
molecule, say a DNA fragment, is pulled by an optical 



trap, the external parameter A corresponds to the trap 
position and W is the work performed on the system by 
the moving trap. Since the force exerted on the molecule 
by the optical trap can be measured for every trap posi- 
tion, the work W can be calculated from the experimental 
data and the average from Equ. ([T]) can be determined 
by repeating the pulling experiment many times. 

While such an analysis of the experimental data is pos- 
sible, one is rarely interested in the free energy of the en- 
tire system (molecule plus trap) as a function of the trap 
position. Rather, the interest usually lies in the equilib- 
rium free energy of the untrapped molecule as a func- 
tion of some conformational degree of freedom such as 
the end-to-end distance. A way to extract that informa- 
tion from non-equilibrium work data has been recently 
proposed by Hummer and Szabo 0, @|. This method 
is based on the insight that the equilibrium distribution 
from which the transformation is initiated can be recon- 
structed from the non-equilibrium distribution generated 
in the process. This reconstruction, in which the bias 
of the trap is removed, requires only knowledge of the 
work carried out in the transformation as well as the 
time-dependent trap potential, both quantities that are 
accessible experimentally. 

Information extracted from single molecule experi- 
ments can be complemented with corresponding non- 
equilibrium computer simulations which provide detailed 
atomistic pictures of the molecule's response to the me- 
chanical perturbation 0, 0, Hcj |. As is the case for the 
experiments, the calculation of equilibrium free energies 
from such simulations using Jarzynski's method or Hum- 
mer and Szabo's refinement becomes increasingly diffi- 
cult if the external perturbation strongly removes the 
system from equilibrium. In this regime, the exponential 
average from Equ. (|TJ) is dominated by a few rare tra- 
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jectories leading to large statistical uncertainties in the 
free energy estimate [ll|. This difficulty may be over- 
come in computer simulations by generating trajectories 
yielding rare but important work values with enhanced 
likelihood, and several techniques to do that have been 
devised [H, 0, 0, EH- Using these biased sampling 
methods, the number of non-equilibrium trajectories re- 
quired to calculate free energies with a given accuracy 
can be dramatically lowered. 

A different type of efficiency enhancement of such non- 
equilibrium fast-switching simulations can be achieved 
by reducing the computational cost required to generate 
the individual trajectories. We have recently proposed a 
method to do that based on the integration of the equa- 
tions of motion with unusually large time steps (l6l . Il7| . 
Since the number of steps required to propagate the sys- 
tem for a given time decreases with increasing time step, 
this approach promises computational savings. Although 
approximate trajectories obtained with large time steps 
mimic the true dynamics of the system only crudely, the 
resulting free energies are in principle exact, the obvious 
limitation being the numerical stability of the integrator. 
In this paper, we investigate how this large-time step 
formalism performs in the simulation of single molecule 
pulling experiments. To do that we first show that the 
Hummer-Szabo procedure can be applied without major 
changes to trajectories obtained with large time steps. 
We then analyse the efficiency of the method for the force 
induced unfolding of deca-alanine and calculate the free 
energy as a function of the end-to-end distance. These 
simulations indicate that for time steps ranging from a 
conservative time step of 0.5 fs to a time step of 3.2 fs, 
just short of the stability limit, the computational ef- 
ficiency is essentially constant. The time steps typically 
used in simulations of such systems are in the range of 1— 
2 fs. A detailed analysis of the trajectories obtained with 
different time steps reveals clear indications of their ap- 
proximate character. In particular, we observe a sudden 
temperature drop of up to 5% during the first few steps 
of the pulling trajectories, as the time step is increased. 
This effect is most likely because of the so-called shadow 
Hamiltonian that is conserved for symplectic integrators 
such as the velocity- Verlet algorithm used in this study. 

The remainder of the paper is organised as follows. 
First, in Sec. [TTJ, we review our large timestep fast- 
switching method and show how the Hummer-Szabo pro- 
cedure can be applied to this case. After that we describe 
the model system in Sec. IIIII Results are presented in 
Sec. IIVI followed by conclusions in Sec. fVl 



II. LARGE TIME STEP FAST SWITCHING 

Any deterministic time evolution in phase space can 
be viewed as a mapping that takes the initial point of 
a trajectory into its final point. Expanding on this per- 
spective, we have recently derived a generalised version 
of the Jarzynski identity that permits to calculate exact 



free energies from approximate large time step trajecto- 
ries both for deterministic as well as stochastic dynamics 
fifty . The basis for this approach is provided by the fol- 
lowing identity that can be shown to hold for a system 
with parameter dependent Hamiltonian 7i(x, A) and a 
general invertible and differentiable mapping <p(x) acting 
on the phase space point x: 
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Here, the angular brackets denote an average over the 
equilibrium (canonical) distribution of the initial state. 
The work function W$, defined as 



Wt(x) = H(<j>(x),X B ) - H(x, A A ) - /T 1 In 



d<j)(x) 



Ox 



(3) 



includes the total energy change caused by switching the 
parameter A from its initial value Xa to its final value 
Xb- In addition, a term depending on the Jacobian 
\d(j)(x) / dx\ of the mapping contributes to the work func- 
tion. If the mapping <fi{x) conserves the canonical distri- 
bution, this entropic term, which takes into account the 
expansion or compression of pha se space volume, corre- 
sponds to the heat exchange In this case, Equ. §3§ 
is equivalent to the first law of thermodynamics. 

To apply Equ. ([3]) to the case of large time step tra- 
jectories, one has to consider the mapping that results 
from a concatenation of molecular dynamics steps. The 
work function then consists of the difference in total en- 
ergy between the initial and final point of the trajectory 
and a sum of terms depending on the phase space vol- 
ume change during each time step. We emphasise that 
in principle Equ. ([3]) yields exact free energy differences 
regardless of the size of the time step. 

The application of the large time step method be- 
comes particularly simple if the algorithm used for the 
integration of the equations of motion conserves phase 
space volume. In this case, the Jacobian is strictly unity 
for each time step and the work function reduces to 
the energy difference accumulated along the trajectory, 
W^(x) = H((j)(x), Ab) — H(x,Xa)- In the present work 
we will only use the velocity- Verlet algorithm, which 
is volume-preserving [T|| [2(j. A detailed discussion of 
the large time step formalism for non volume-preserving 
integrators and related complications can be found in 
Ref. P3. 

Since Jarzynski's identity holds for a time evolution 
resulting from the integration of the equations of motion 
with large time step, one expects that also the Hummer- 
Szabo procedure mentioned above remains valid in this 
case. We next show that this is indeed the case. Consider 
a system with Hamiltonian 



H(x,X B )=H (x) + V(q{x);X), 



(4) 



where TLq{x) is the molecular Hamiltonian and V(q(x), A) 
is the parameter dependent trap potential coupling to the 
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variable q(x). What one would like to calculate is the free 
energy F(q) defined up to a constant as 



-™ = (<%-<z(z)]> , 



(5) 



where the brackets (•••)o denote an average over the 
equilibrium distribution of the molecular system without 
the trap, i.e., over the distribution 
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The goal now is to derive an expression that permits 
to compute the free energy from the phase space distri- 
bution resulting from application of the mapping <f>(x) 
where at the same time the parameter A is changed from 
its initial value Xa to its final value A#. To do so, we 
first perform a transformation of variables from x to 
y = <j)~ l (x) obtaining 



dy 
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where Qq = J dx exp[— /3Hq(x)] is the partition function 
of the system without trap. Multiplication and division 
of the integrand with exp{— f3H(y, Xa)] ~ (3V(q(y), As)} 
then yields 



where 



Qo 



(8) 
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is the partition function of the entire system, including 
the trap at position Xa and the average is over initial con- 
ditions canonically distributed with respect to 7i(x, Xa)- 
Equation ([5]) implies that the entire free energy profile 
F(q) = — kBThx(6[q — q(x)])o can be calculated up to a 
constant by histogramming the variable q at the end of 
the transformation </>(x) and weighing each contribution 
to the histogram with the work exponential exp(— (iW$). 
This weight essentially takes into account the different 
phase space probability of a particular point in the equi- 
librium ensemble and the ensemble generated by the 
mapping. Equation ((5)1 is analogous to the central re- 
sult of Hummer and Szabo (Equ. (7) of Ref. Q), which 
therefore is valid also for general phase space mappings. 

In the case of a molecular dynamics simulation, the 
mapping <fi(x) consists of a concatenation of a certain 
number of time steps. Equation (jHJ) can be applied at 
each stage of the mapping and, accordingly, we rewrite 
it as 
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Here, the index i refers to the number of steps the map- 
ping <pi and W r j >i consists of is the work accumulated in 



the first i steps. The parameter A is changed in discrete 
steps Xi such that Ao = Xa and A„ = A^ and n is the 
total number of steps. 

In principle, the full free energy profile can be obtained 
from data at one single time. For better accuracy one 
can combine the free energy profiles for all time steps 
using the weighted histogram technique introduced by 
Ferrenberg and Swendsen plj : 
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Combining the histograms for different times in this way, 
each histogram contributes most where it has the high- 
est accuracy. In conclusion, we have demonstrated that 
the method of Hummer and Szabo can be transferred es- 
sentially without changes to the case of large time step 
trajectories. 



III. MODEL AND PULLING PROCESS 

To investigate the computational efficiency of pulling 
simulations with large time steps we study the force- 
induced unfolding of a deca-alanine molecule in vacuum. 
This oligo-peptide with an acetylated N-terminus and an 
amidated C-terminus consists of ten alanine residues and 
forms an a-helix that is stable in vacuum at room tem- 
perature (Fig. [J). It has N = 109 atoms and is therefore 
small enough to permit the calculation of a large number 
of trajectories in a reasonable amount of time, but is, on 
the other hand, close to the systems treated in molecular 
biology or chemistry. Except the differing termini, the 
model is identical to the one studied by Park et. al. [12] • 
By using this system studied previously as our test case, 
we can check the results obtained with large time steps 
not only against our own calculations with "safe" time 
steps, but against an external reference. 





FIG. 1: (Color online) Deca alanine in its initial a-helical 
(top) and its final coil configuration (bottom). The green 
atom on the right hand side stays fixed throughout the simu- 
lation while the green atom on the left hand side moves in a 
harmonic trap that is translated from right to left - as indi- 
cated by the spring. 



4 



The unfolding process from the helix to the coil config- 
uration is induced by imposing a time-dependent trap po- 
tential V(q; t) on the amide nitrogen atom of the capped 
C-terminus at position q while keeping the nitrogen atom 
of the first residue fixed at the origin. The harmonic trap 
potential V(q; t) is given by 



where 



V{q;t) = -[q-z{t)f 



z(i) = Zi + vt. 



(12) 



(13) 



is the trap position and k is the force constant of the 
potential which is set to k = 14.38 kcal/mol A in all 
simulations. The trap is moved from its initial (t = 0) 
position Zi to its final (t = r) position Zf at constant 
speed v = (z f — Zi) / V '. The initial trap position of Zi = 13 
A is chosen such that the system starts in a slightly com- 
pressed state where the end-to-end distance is smaller 
than the equilibrium distance, which lies at about 15.4 
A. The final position Zf = 33 A leads to a fully unfolded 
state. The trap speed v must be chosen sufficiently low 
such that it allows good sampling using the Jarzynski 
equality. On the other hand, it should be high enough 
to permit sampling of sufficiently many trajectories. Fol- 
lowing Park et al. HH, we chose a trap speed of 0.01 
A/ps, which amounts to a trajectory length of 2 ns. 

The canonically distributed initial conditions, as re- 
quired by Jarzynski's theorem, are generated from an 
equilibrium Langevin molecular dynamics simulation at 
z(t) = Zi with a small time step of 1 fs. From this 
"base trajectory" we start a non-equilibrium pulling tra- 
jectory every one thousand time steps using the velocity 
Verlet algorithm. All our simulations were performed 
using the CHARMM force field version 27 with both 
CHARMM and the NAMD simulation package, respec- 
tively mni n hi. 



IV. RESULTS 

In this section we first show that we obtain identical 
free energy profiles (well within statistical error bars) , re- 
gardless of the time step used. We then discuss an inter- 
esting effect we observed for the temperature behaviour 
in our large timestep simulations. Finally, we address the 
central point of this paper, the efficiency of the large time 
step method for the calculation of free energy profiles. 



1. Free energy profiles 

From repeated pulling simulations one can collect sep- 
arate distributions of the end-to-end distance q for every 
time U. Examples of such distributions obtained for the 
final position of the pulling trap for different time steps 
At are depicted in Fig. [2] along with the trap potential. 
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FIG. 2: (Color online) Histograms P(q) of the end-to-end- 
distance q for various timesteps obtained from 5000 trajec- 
tories at t = 2 ns corresponding to the final position of the 
trap. 



Note that the histograms do not differ significantly for 
the different time steps and that the distributions are 
not centred at the minimum of the trap potential but to 
the left of it, indicating that the molecule slightly lags 
behind the trap potential. 

Using the procedure of Hummer and Szabo, we have 
combined the distributions corresponding to different 
times obtaining the free energy as a function of the end- 
to-end distance q. For the application of Equ. (fTTj) it 
is important that the histograms of the end-to-end dis- 
tance have sufficiently narrow bins. If not, the correction 
for the pulling potential in the denominator of Equ. (jlip 
varies too much within a bin, thus giving the wrong cor- 
rection factor. In our simulations we used a bin size of 
0.02 A. The free energy profiles calculated with Equ. (fTTj) 
are depicted in Fig. [3] and agree well with the curves cal- 
culated by Park et. al. [22|. As can be inferred from the 
figure, all time steps yield the same result within statis- 
tical deviations. The free energy profile F(q) displays a 
minimum at an end-to-end distance of q ss 15 A corre- 
sponding to the a-helical native state. For increasing q, 
the free energy then grows almost linearly up to q ~ 25 A. 
In this regime, the hydrogen bonds in the helix are con- 
secutively cleaved. After all hydrogen bonds have been 
broken, the free energy keeps growing, but at a smaller 
rate. No stable state exists for the unfolded molecule. 



2. Artifacts of large time step trajectories 



The instantaneous temperature 

2 (l?kin) 
3Nk B 



Ti = 



(14) 



during the pulling process averaged over two-thousand 
trajectories is depicted in Fig. 3] for various time steps. 
The overall temperature behaviour is similar for all time 
steps: after a small rise occurring at a trap position of 
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FIG. 3: (Color online) Free energy profiles F(q) as function 
of the end-to-end distance q for various timesteps obtained 
from 5000 trajectories using the Hummer-Szabo procedure. 
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FIG. 4: (Color online) Instantaneous temperature T, as a 
function of the trap position z obtained for various timesteps 
and averaged over 2000 trajectories. 



about 15 A, the temperature decreases steadily through- 
out the switching process. However, the temperature 
traces differ by an essentially constant offset which is due 
to a very sharp decline of up to five percent from the 300 
K of the Langevin base trajectory. This initial tempera- 
ture drop, which occurs during the first few integration 
steps, grows with increasing time step and is accompa- 
nied by an increase in the potential energy such that the 
total energy remains approximately constant. The po- 
tential energy increase is mainly caused by an increase 
of the angle bending potential energy terms and, to a 
smaller extent, of the bond stretching interactions. 

The redistribution of energy from kinetic to potential 
terms suggests that the equipartition theorem may be 
violated for large time step trajectories. We indeed find 
that in straightforward equilibrium MD simulations car- 
ried out with time steps larger than 2 fs the kinetic en- 
ergy of the hydrogen atoms is up to w 10 percent less 
than that of heavier atoms. 

This phenomenon is most likely linked to the nearly 
harmonic potentials used to model the angle bending 
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FIG. 5: Schematic phase space portrait of the exact dynam- 
ics of a harmonic oscillator (solid line) and of the dynamics 
obtained with large time steps (dashed line). 



and bond stretching interactions and can be made plausi- 
ble by considering a one-dimensional harmonic oscillator 
with Hamiltonian 



p 2 moj 2 q 2 
U{q,p) = — + — - — . 
2m 2 



(15) 



Here, q and p denote the position and the velocity of 
the oscillator and m and u> are its mass and frequency, 
respectively. As is well known [13 , |28| . l29l |30I | , symplec- 
tic integration schemes like the velocity Verlet algorithm 
do not conserve the actual Hamiltonian of a system, but 
a timestep-dependent "shadow Hamiltonian" Tip which 
converges to the real Hamiltonian only in the limit of 
At — > 0. For the harmonic oscillator this shadow Hamil- 
tonian is known analytically [3lj ] : 
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Starting from a given phase space point (qo,po) the nat- 
ural dynamics of the harmonic oscillator traces out an 
ellipsoid on which the Hamiltonian from Equ. (|15[) is con- 
served (solid line in Fig. [5]). However, a trajectory ob- 
tained by integrating the equations of motion of the har- 
monic oscillator with finite time step At lies on a curve 
on which the shadow Hamiltonian TLp(x, v) rather than 
the Hamiltonian TC is conserved (dashed line in Fig. [5]) . It 
follows from the specific form of the shadow Hamiltonian, 
that this curve, which is also an ellipsoid, is compressed 
on the momentum axis by a factor 



pi + m 2 LU 2 q 2 (l - w 2 Ai 2 /4) 
p\ + m 2 uj 2 ql 

but expanded on the position axis by a factor 

p 2 (l - w 2 At 2 /A)- 1 + m 2 tu 2 q 2 
Pq + m 2 uj 2 q 2 



(17) 



(18) 



with respect to the ellipsoid defined by a constant Hamil- 
tonian. Thus, increasing the integration time step At 
leads to a decrease in the average kinetic energy and, at 
the same time, to an increase of the average potential 
energy. 

Under the assumption that the dynamics obtained by 
application of the velocity Verlet algorithm is nearly con- 
tinuous, a canonical average of the kinetic and potential 



6 



energy averaged over one oscillation period yields 



and 
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u 2 At 2 



k B T f l-uj 2 At 2 /8 \ 
l-uj 2 At 2 /A) 



(19) 
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respectively. The deviation AT^i n /T of the kinetic tem- 
perature Tkin = 2(K)/kB from the temperature T of the 
distribution of the initial conditions is then obtained from 
the averaged kinetic energy according to: 



ATki 
T 



n 



T 



uo 2 At 2 



T 



(21) 
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FIG. 6: Relative deviation AT/T of the kinetic temperature 
from the temperature of the initial conditions (symbols) as a 
function of At 2 . The line is a linear fit to the data with slope 
-0.0043 fs~ 2 . 

The quadratic dependence of the temperature drop 
from the time step expected for a system with nearly 
harmonic degrees of freedom can be tested for the deca- 
alaninc molecule. Figure [6] shows the relative tempera- 
ture drop obtained from the data depicted in Fig. 0] as 
a function of At 2 along with a linear fit to the data. 
Approximately, the relative temperature drops follow 
the quadratic behaviour found in the harmonic oscilla- 
tor. Moreover, the slope of the fitted line is consis- 
tent with the assumption that each hydrogen atom is 
involved in exactly one harmonic degree of freedom and 
the frequency of the corresponding motion is of the order 
v w 3000cm- 1 . 

The above considerations suggest that the sudden drop 
of the kinetic temperature observed for the deca-alanine 
with increasing time step may be related to the particular 
form of the shadow Hamiltonian corresponding to the 
nearly harmonic degrees of freedom of the molecule. We 
stress, however, that this effect has no impact on the 
large timestep method, because the only requirements are 
a canonical ensemble of initial conditions and a volume 
preserving integrator such as the velocity Verlet scheme. 



As long as these conditions are met, the large time step 
method is valid regardless of the unusual behaviour of 
some energy components. 



3. Efficiency 

To assess the efficiency of the large timestep algorithm 
we need a way to determine the statistical error in the free 
energy profile as a function of the number of trajectories 
used in the Jarzynski average. For this purpose we use 
a block av erag ing method [32j, analogous to Zuckerman 
and Woolf [27]. First we partition all N trajectories into 
M "blocks" of size n <C N and introduce the free energy 
F^(q) calculated by applying Equ. (fTTj) only to the n 
trajectories of block i. Then, wc compute the statistical 
deviation - Equ. (f2"2"| - of F^(q) from the mean F(q), 
averaged over all trajectories. To improve the accuracy 
of this error estimate we average also over the end-to-end 
distances q: 



C 2 (n) = 



(22) 



Here, S is the number of end-to-end distances over which 
the average extends and in all our calculations S = 1000. 
The quantity C 2 (n) is the mean square deviation of the 
block free energy from the free energy obtained from all 
trajectories. For a different approach to efficiency es- 
timation in fast switching simulations, which, however, 
requires a larger computational effort, see [lq ]. 

From the resulting curves C 2 (n), which are depicted 
in Fig. [7J we can extract the number N&t of trajectories 
needed to obtain an accuracy of k^T. To do that we 
determine the intersection of the correlation curves C 2 (n) 
with the error threshold, shown as a red dashed line in 
Fig. [7J The threshold of k&T is arbitrary, but, as the 
inset of Fig. [JJ indicates, the relative numbers N^t for 
different timesteps are largely independent of the choice 
of the threshold. 

The total computational cost required to obtain a free 
energy profile with an accuracy of k^T is proportional to 
the number of force evaluations performed in the calcu- 
lation and thus to the number L of time steps per tra- 
jectory: 



cost (At) = LN At 



At' 



-N 



Al- 



(23) 



The cost function cost(At) is the total number of molec- 
ular dynamics steps required to calculate the free energy 
with accuracy k^T, or, in other words, the total com- 
putational cost of the simulation in units of the cost of 
one single molecular dynamics step. The behaviour of 
the computational cost as a function of the integration 
time step At (see Fig. [5]), indicates that in the range 
At = 1.0 — 3.2 fs the computational cost of the free en- 
ergy calculation is essentially constant while time steps 
shorter than « 1 fs lead to an increased computational 
cost. 
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FIG. 7: (Color online) C 2 (n) as a function of the block 
size n for different timesteps and the threshold (ksT) 2 ~ 
0.355(kcal/mol) 2 (red dashed line). The inset shows the same 
graphs on a doubly logarithmic scale. 
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calculation of valid free energies from computationally 
inexpensive trajectories. In the present work we have 
adapted the Hummer-Szabo procedure to large time step 
dynamics and have shown that it is applicable to the 
simulation of complex biomolecules. 

Approximate large time step trajectories generated 
with the velocity Verlet algorithm display an interesting 
initial temperature drop, which is most likely due to the 
specific properties of the shadow Hamiltonian conserved 
by the discrete dynamics. Although this effect leads to 
inaccurate trajectories violating the equipartition theo- 
rem, it is irrelevant for the validity of the large timestep 
formalism since all its requirements, a phase space vol- 
ume conserving integrator and a canonical distribution 
of initial conditions, are met also in this case. It should 
be noted that the temperature jump indicating the ap- 
proximative nature of the trajectories can be observed 
already at a time step of 2 fs, a value which is often 
used in equilibrium simulations of protein and peptide 
systems. Thus, one may argue that the large time step 
approach is more correct than a comparable calculation 
based on equilibrium methods alone. 

The potential computational benefit of the large time 
step method is in part compensated by a growth of sta- 
tistical fluctuations observed in simulations using compu- 
tationally inexpensive trajectories generated with larger 
time steps. All time steps from about I fs up to the sta- 
bility limit of the integrator yield essentially the same 
efficiency Time steps smaller than about f fs, however, 
lead to an efficiency loss. 



FIG. 8: Computational cost as a function of time step At. 
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